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PACS 45 . 70 . -n - Granular systems 

PACS 02 . 50 . Ey - Stochastic processes 

PACS 05 . 40 . -a - Fluctuation phenomena, random processes, noise, and Brownian motion 

Abstract. - A Generalized Langevin Equation with exponential memory is proposed for the 
dynamics of a massive intruder in a dense granular fluid. The model reproduces numerical cor- 
relation and response functions, violating the equilibrium Fluctuation Dissipation relations. The 
source of memory is identified in the coupling of the tracer velocity V with a spontaneous local ve- 
locity field U in the surrounding fluid: fluctuations of this field introduce a new timescale with its 
associated lengthscale. Such identification allows us to measure the intruder's fluctuating entropy 
production as a function of V and U , obtaining a neat verification of the Fluctuation Relation. 
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Models of granular fluids are a natural framework where 
the issues of non-equilibrium statistical mechanics can 
be addressed pQ. Due to dissipative interactions among 
the microscopic constituents, energy is not conserved and 
external sources are necessary in order to maintain a 
stationary state. Heat fluxes and currents continuously 
ass through the system, time reversal invariance is bro- 
ken and consequently, properties such as the Equilibrium 
Fluctuation-Dissipation relation (EFDR) do not hold. In 
recent years, a rather complete theory, at least in the di- 
lute limit, has been developed and numerous aspects have 
been clarified, in good agreement with numerical simula- 
tions 0E1- However, a general understanding of dense 
granular fluids is still lacking. A common approach is the 
so-called Enskog correction [2B] , which reduces the break- 
down of Molecular Chaos to a renormalization of the col- 
lision frequency. In cooling regimes, the Enskog theory 
may describe strong non-equilibrium effects, due to the 
explicit cooling time-dependence [5]- However it cannot 
describe dynamical effects in stationary regimes, such as 
large violations of the Einstein relation [HIE] • 



In this letter, we propose a model for the dynamics of 
a massive tracer moving in a gas of smaller granular par- 
ticles, both coupled to an external bath. In particular, 
taking as reference point the dilute limit, where the sys- 
tem has a closed analytical description [8], we suggest a 
Generalized Langevin Equation (GLE) with an exponen- 
tial memory kernel as first approximation capable of de- 
scribing the dense case. Here, the main features are: i) the 
decay of correlation and response functions is not simply 



exponential and shows backscattering [51[TU] and ii) the 
EFDR [TT1II2 of the first and second kind do not hold. 
In the model we propose, detailed balance is not necessar- 
ily satisfied, non-equilibrium effects can be taken into ac- 
count and the correct behavior of correlation and response 
functions is predicted. Furthermore, the model has a re- 
markable property: it can be mapped onto a two-variable 
Markovian system, i.e. two coupled Langevin equations 
with simple white noises. The auxiliary variable can be 
identified in the local velocity field spontaneously appear- 
ing in the surrounding fluid. This allows us to measure 
the fluctuating entropy production |13j . and fairly verify 
the Fluctuation Relation [T^HHHS] . This is a remarkable 
result, if considered the interest of the community |16| and 
compared with unsuccessful past attempts [TTHT5] . 

We consider an "intruder" disc of mass mo = M and 
radius R, moving in a gas of N granular discs with mass 
rrii = m (i > 0) and radius r, in a two dimensional box of 
area A — L? . We denote by n — N/A the number den- 
sity of the gas and by <j) the occupied volume fraction, i.e. 
<f> = ir(Nr 2 + R 2 )/A and we denote by V (or Vq) and v 
(or Vi with i > 0) the velocity vector of the tracer and 
of the gas particles, respectively. Interactions among the 
particles are hard-core binary instantaneous inelastic col- 
lisions, such that particle i, after a collision with particle 
j, comes out with a velocity 



1 



■[(Vi -Vj) ■h]fi 



(1) 



where n is the unit vector joining the particles' centers of 
mass and a £ [0, 1] is the restitution coefficient (a = 1 



p-1 



A. Sarracino, D. Villamaina, G. Gradenigo A. Puglisi 



is the elastic case). The mean free path of the intruder 
is proportional to fo = l/(n(r + R)) and we denote by 
t c its mean collision time. Two kinetic temperatures 
can be introduced for the two species: the gas granular 
temperature T g = m(v 2 )/2 and the tracer temperature 
T tr = M(V 2 )/2. 

In order to maintain a granular medium in a fluidized 
state, an external energy source is coupled to each particle 
in the form of a thermal bath [T9H2T] (from hereafter, ex- 
ploiting isotropy, we consider only one component of the 
velocities) : 



miVi(t) = -7&Vi(t) + fi(t) + £ b (t). 



(2) 



Here /,(£) is the force taking into account the collisions 
of particle i with other particles, and £&(£) is a white 
noise (different for all particles), with (£, b {t)} = and 
(&(*)&(*0) = 1T blb 5{t - t'). The effect of the external 
energy source balances the energy lost in the collisions 
and a stationary state is attained with rrii(v 2 ) < T b . 

For low packing fractions, <fi < 0.1, and in the large 
mass limit, m/M <C 1, using the Enskog approximation 
it has been shown [8] that the dynamics of the intruder is 
described by a linear Langcvin equation. In this limit the 
velocity autocorrelation function shows a simple exponen- 
tial decay, with characteristic time M/Te, where 



lb+ln 



with 



1 E 



92{r + R) 



^27rmT g (l + a) 

(3) 

and gi(r + R) is the pair correlation function for a gas 
particle and the intruder at contact. Time- reversal and 
the EFDR, which are very weakly modified for uniform 
dilute granular gases [5J[221[23] , become perfectly satisfied 
for a massive intruder. The temperature of the tracer is 
computed as T$ = (~y b T b + ^^T g )/T E . For a general 
study of a Langevin equation with "two temperatures" 
but a single time scale (which is always at equilibrium), 
see also [23] . 

As the packing fraction is increased, the Enskog approx- 
imation is less and less effective in predicting the memory 
effects and the dynamical properties of the system. In par- 
ticular, velocity autocorrelation C( t) = (V (t)V (0)) / (V 2 ) 
and linear response function R(t) = 5V(t)/SV(0) (i.e. the 
mean response at time t to an impulsive perturbation ap- 
plied at time 0) present an exponential decay modulated 
in amplitude by oscillating functions [TU]. Moreover vio- 
lations of the EFDR C(t) = R(t) (Einstein relation) are 
observed for a < 1 [71125]. 

Molecular dynamics simulations of the system have been 
performed by means of a standard event driven algorithm 
to treat hard core interactions: the algorithm is supple- 
mented with a "driving event" at times which are multi- 
ples of a small timestep (smaller than all timescales) which 
update the velocity of all particles by a discretized version 
of Eq. ([2]). In the simulations we have measured C(t) and 
R(t), for several different values of the parameters a and 
</>. In Fig. [TJ symbols correspond to the velocity correla- 




Fig. 1: (Color online). Semi-log plot of C(t) (symbols) for 
different values of 4> — 0.01, 0.1, 0.2, 0.33 at a = 0.6. Times are 
rescaled by the mean collision time r c . Continuous lines are 
the best fits obtained with Eq. (J9]). Inset: C(t) and the best 
fit in linear scale for d> = 0.33 and a — 0.6. 



tion functions measured in the inelastic case, a = 0.6, for 
different values of the packing fraction <f>. The other pa- 
rameters are fixed: N = 2500, m = 1, M = 25, r = 0.005, 
R = 0.025, T b = 1, 7 6 = 200. 

Notice that the Enskog approximation [211] cannot pre- 
dict the observed functional forms, because it only modi- 
fies by a constant factor the collision frequency. In order to 
describe the full phenomenology, a model with more than 
one characteristic time is needed. As a first proposal, we 
consider a Langevin equation with a single exponential 
memory kernel [26,27 



MV(t) = dt' T(t - t')V{t') + £'(£), 



where 

r(£) = 2 7o £(t)+7i/Tie- t / Tl 

and £'(t) = £ (t) + £i(t), with 

(£ (t)£ (t')) = 2T o7o 5(t-t'), 
(£ 1 (t)£ 1 (t'))=T l7l /r 1 e-^- t '^ 



(4) 
(5) 

(6) 

(7) 



and (£i(i)£o(t')) = 0. In the limit a -> 1, the parameter 
Xi is meant to tend to T in order to fulfill the EFDR 
of the 2nd kind (£' {t)£' (?)) = T T(t - £')■ Within this 
model the dilute case is recovered if 71 — >• 0. In this limit, 
the parameters 70 and To coincide with Te and T® of the 
Enskog theory [8]. 

The exponential form of the memory kernel can be jus- 
tified within the mode-coupling approximation scheme. In 
this framework [28) , it can be written as a sum of two con- 
tributions: r(t-i') = p 1 8{t-t')+p 2 ^ (t-f), where ft and 
P2 are model dependent coefficients, and T(t — t') is a sum 
over modes q oip(q)e~^ +D ^ q (*~* ', where p{q) weights the 
modes relevant for the dynamics of the tracer. Here D and 
v are the diffusion coefficient and the kinematic viscosity of 
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Table 1: Parameters of model (|10|1 . as obtained by fitting the numerical data (see text for details). 
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the fluid respectively. Following an old recipe [26], tested 
with success in equilibrium contexts, we assume that, for 
not too high packing fractions, memory arises due to re- 
collisions within a limited region at distance ~ Ai around 
the tracer and that this can be modeled by an effective 
p(q) which is peaked around q\ — 2ir/\i, i.e. a single mode 
contributes to the sum, yielding f (t — t') ~ e~ ( - l,+D ^ qi( - t ~ t ' s> 
and then 



ri = \\{2tt)~ 2 {v + D) 



^(Ai/O 



<J\2 



(8) 



t,? and Iq the fluid mean free time and mean free 



with 

path respectively. Eq. (O relates the time-scale n, char- 
acterizing the tail of the memory kernel, with a typical 
length-scale Ai present in the system. This length-scale 
will turn out to play a central role in the following. 
The model 0} predicts C = f c (t) and R = f R (t) with 



fc(R) = e gt [cos(ujt) + ac{R) sin(wt)]. 



(9) 



g, w, ac and clr are known algebraic functions of 70, 
Tq, 71, T\ and T±. In particular, the ratio ac/a>R = 
[T - n(Ti - To)]/[T + n(Ti - T )], with fi = 7l /[( 7o + 
7i)(7o/Mn— 1)]. Hence, in the elastic (T. — > To) as well as 
in the dilute limit (71 — > 0), one gets ac = cir and recov- 
ers the EFDR C(t) — R(t). In Fig. [T] the continuous lines 
show the result of the best fits obtained using Eq. ([9]) for 
the correlation function, at restitution coefficient a = 0.6 
and for different values of the packing fraction cf>. The 
functional form fits very well the numerical data. 

Looking for an insight of the relevant physical mecha- 
nisms underlying such a phenomenology and in order to 
make clear the meaning of the parameters, it is useful to 
map Eq. ([!]) onto a Markovian equivalent model by intro- 
ducing an auxiliary field |29j : 



MV = 

(J = 



-lo{V -U) + ^2T Q -fo£v 

---—v + xl^z 

n 70 n V 7o r i 



(10) 



where £ y and £ u are white noises of unitary variance. The 
variable 

U(t) oc 71/(7170) / e-^fV(i') + £i(i')]<ft' (11) 



is determined up to a multiplicative factor, as it can be 
checked by direct substitution. In the chosen form <jTUJ) , 
the dynamics of the tracer is remarkably simple: indeed V 
follows a memoryless Langcvin equation in a Lagrangian 
frame with respect to a local field U. In the dilute limit 
this is exact (see Appendix of [8]) if C/ is the local average 
velocity field of the gas particles colliding with the tracer. 
Extrapolating such an identification to higher densities, 
we are able to both assign a meaning and predict a value 
for most of the parameters of the model: 1) the self drag 
coefficient of the intruder in principle is not affected by 
the change of reference to the Lagrangian frame, so that 
7o ~ T-e; 2) for the same reason To ~ T tr is roughly the 
temperature of the tracer; 3) 7*1 is the main relaxation 
tinrc of the average velocity field U around the Brown- 
ian particle; 4) 71 is the intensity of coupling felt by the 
surrounding particles after collisions with the intruder; 5) 
finally T\ is the "temperature" of the local field U, eas- 
ily identified with the bath temperature T\ ~ T,: indeed, 
thanks to momentum conservation, inelasticity does not 
affect the average velocity of a group of particles almost 
only colliding among themselves. 

To find a confirmation of the above hypothesis, we have 
explored the region of the space of parameters a £ [0.6, 1] 
and £ [0.01,0.33]. From the simultaneous fit of the 
numerical data for correlation and response functions 
against Eqs. (O we can determine the set of parameters 
{g , oj , ac 1 o-r, (V 2 )}- Then, by inverting the relations be- 
tween them and the set {70, To, 71, ri,Ti}, we are even- 
tually able to determine all the parameters entering Qj. 
In Table Q] such values are reported, together with the 
predictions given by the Enskog approximation (last four 
columns). The statistical error on these values is about 
1%. We used the external parameters mentioned before, 
changing a or the box area A (to change </>): this makes 
the limit — > equivalent to j g ~ l/l — > ("super- 
dilute" limit). The last row reports about the true dilute 
limit: i.e. R is reduced, at fixed Iq (equal to the value of 
the previous case <j> = 0.2), in order to get = 0.01 and 
7 9 > 0. Notice that in the two dilute cases the simple 
Langevin equation is recovered (71 = 0) and the depen- 
dence on the parameter T\ disappears. Remarkably our 
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predictions 70 ~ r^, To ~ T tr and Ti ~ T, are fairly 
verified. The coupling time rj increases with the pack- 
ing fraction and, weakly, with the inelasticity. In the most 
dense cases it appears that 71 ~ 7,? ex <f>: this is confirmed 
in the "super-dilute" limit, but cannot hold in the dilute 
one, where 71 — >• <C 7J 5 . It is also interesting to notice 
that at high density T tr ~ T g ~ T^ , which is probably due 
to the stronger correlations among particles. Finally we 
notice that, at large (j>, T tr > T t f, which is coherent with 
the idea that correlated collisions dissipate less energy. 

A fundamental feature of this model is its ability to 
reproduce violations of EFDR. In Fig. [5J we plot the cor- 
relation and response functions in a dense case (elastic and 
inelastic): symbols correspond to numerical data and con- 
tinuous lines to the best fit curves. In the inelastic case, 
deviations from EFDR R(t) = C(t) are clearly observed. 
In the inset of Fig. |5] the ratio R(t)/C(t) is also reported. 
It is interesting to note that a relation between the re- 
sponse and correlations measured in the unperturbed sys- 
tem still exists, but - in the non-equilibrium case - must 
take into account the contribution of the cross correlation 
(V(t)U(0)), i.e.: 



R(t) = aC(t) + b(V(t)U(Q)) 



(12) 



with a = [1 - 7i/M(T - Tx)Ct a ] and b = (T - Ti)ft b , 
where H a and Of, are known functions of the parameters 
(see for instance [21] )• At equilibrium, where To = T\, the 
EFDR is recorvered. 

The mathematical definition of the auxiliary variable 
U, Eq. (fTTj) . which requires the knowledge of a part of the 
noise £\, makes it very difficult to be measured in simu- 
lations or in experiments. But the above discussion has 
shown that U represents a spontaneous local velocity field 
interacting with the tracer: therefore it can be measured 
in the following manner. We fix a distance I and average 
the velocity of the gas particles within a circle Ci of ra- 
dius I + R centered on the tracer. In this way we define 
Ui = 1/A^^iec. Vil wnere Ni is the number of particles 
in Ci. Two methods are available to estimate the correct 
length I*, which is difficult to be predicted on a general 
ground. A first guess is provided by identifying it with 
Ai, which can be obtained by inverting Eq. (|8]) after hav- 
ing measured n , using the known values of D and v in a 
granular fluid. The second method is to measure the cor- 
relations (VUi) and (Uf) and find the best value l cor such 
that (VUi cor ) ~ (VU) and (E/^J ~ (U 2 ) (where (VU) 
and (U 2 ) are easily computed from the model, once all the 
parameters have been determined fitting C(t) and R(t)). 
Remarkably, the two estimates give compatible results and 
identify a narrow range of values for I* ~ Ai ~ l cor - Hence, 
one can identify U ~ Ui* and the auxiliary variable can 
be directly measured in numerical simulations and exper- 
iments. 

An important independent assessment of the effective- 
ness of model (jlj) comes from the study of the fluctuat- 
ing entropy production |13j which quantifies the deviation 




Fig. 2: (Color online) . Correlation function C(t) (black circles) 
and response function R(t) (red squares) for a = 1 and a — 
0.6, at <j> — 0.33. Continuous lines show the best fits curves 
obtained with Eqs. (J9]). Inset: the ratio R(t)/C(t) is reported 
in the same cases. 



from detailed balance in a trajectory. Given the trajectory 
in the time interval [0, t], {^(s)} , and its time-reversed 
{IV{s))l = {-V(t - s)} , in Ref. [3D] it has been shown 
that the entropy production for the model (@| takes the 
form 



log 



P({V(s)Y Q ) 



7o 



1 



1 



ds V{s)U(s) 



P({TV(s)Yo) 

(13) 

Boundary terms - in the stationary state - are sublead- 
ing for large t and have been neglected. This functional 
vanishes exactly in the elastic case, a = 1 , where equipar- 
tition holds, T\ = To, and is zero on average in the dilute 
limit, where (VU) = 0. Formula (|T3|) reveals that the lead- 
ing source of entropy production is the energy transferred 
by the "force" jqU on the tracer, weighed by the differ- 
ence between the inverse temperatures of the two "ther- 
mostats" . Following the procedure described above, in the 
case <p = 0.33 and a = 0.6, we estimate for the correla- 
tion length I* ~ 9r ~ 61q. Then, measuring the entropy 
production of Eq. (fT3|) (by replacing U(t) with Ui>) along 
many trajectories of length t, we can compute the proba- 
bility P(St = x) and compare it to P(£< = —x), in order 
to verify the Fluctuation Relation 



log 



P(Zt 



P(£ t = -x) 



(14) 



In Fig.[3Jwe report our numerical results. The main frame 
confirms that at large times the Fluctuation Relation (fT4"]) 
is well verified within the statistical errors. The inset 
shows the collapse of logP(S t )/i onto the large devia- 
tion rate function for large times. Notice also that for- 
mula (|T3| does not contain further parameters but the 
ones already determined by correlation and response mea- 
sure, i.e. the slope of the graph is not adjusted by further 
fits. Indeed a wrong evaluation of the weighing factor 
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Fig. 3: (Color online). Check of the fluctuation relation (|14|) 
in the system with q = 0.6 and <j> = 0.33. Inset: collapse of 
the rescaled probability distributions of £4 at large times onto 
the large deviation function. 



{I /T - 1/Ti) w (1/T tr - 1/T 6 ) or of the "energy injection 
rate" joU(t)V(t) in Eq. (|T5|) could produce a completely 
different slope in Fig. [3J 

In conclusion, we designed a first granular dynami- 
cal theory describing non-equilibrium correlators and re- 
sponses for a massive tracer. The value of this proposal 
is to offer a significant insight into the mechanisms of re- 
collision and dynamical memory and their unexplored re- 
lation with the breakdown of equilibrium properties. It 
is remarkable that velocity correlations (V(t)U(t')) be- 
tween the intruder and the surrounding velocity field are 
responsible for both the violations of the EFDR and the 
appearance of a non-zero entropy production, provided 
that the two fields are at different temperatures. Small 
non-Gaussian corrections [53], always present in granular 
fluids, are neglected here in favor of the largest contribu- 
tion given by memory terms to violations of EFDR and 
entropy production. For some of the parameters in the 
theory (70 ~ Te, To ~ T tr and T± ~ Tf,) we have reason- 
able predictions, while n and 71, related to the coupling 
between U and V, deserve further investigations. Close 
analytical predictions of all the parameters could be ob- 
tained through a full kinetic theory (beyond Enskog), also 
to deduce eventual extensions to the case M ~ to, larger 
densities, and hard spheres. 
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